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ABSTRACT: Graph i cal Gauss i an m odels with edge and vertex symmetries were intro¬ 
duced by IHoisgaard fc Lauritzenl (20081 ] who also gave an algorithm to compute the maxi¬ 
mum likelihood estimate of the precision matrix for such models. In this paper, we take a 
Bayesian approach to the estimation of the precision matrix. We consider only those models 
where the symmetry constraints are imposed on the precision matrix and which thus form 
a natural exponential family with the precision matrix as the canonical parameter. 

We first identify the Diaconis-Ylvisaker conjugate prior for these models and develop a 
scheme to sample from the prior and posterior distributions. We thus obtain estimates of 
the posterior mean of the precision matrix. 

Second, in order to verify the precision of our estimate, we derive the explicit analytic 
expression of the expected value of the precision matrix when the graph underlying our model 
is a tree, a complete graph on three vertices and a decomposable graph on four vertices with 
various symmetries. In those cases, we compare our estimates with the exact value of the 
mean of the prior distribution. We also verify the accuracy of our estimates of the posterior 
mean on simulated data for graphs with up to thirty vertices and various symmetries. 

KEY WORDS: Conditional independence, symmetries, covariance estimation, trees, Diaconis- 
Ylvisaker conjugate priors, Metropolis-Hastings. 
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1 Introduction 


Given an undirected graph G = (V, E ) where V is the set of vertices and E C V x V 
is the set of undirected edges denoted (i. j), a graphical Gaussian model is a family of 
Gaussian distribution for X = (X v , v £ V) where the conditional independences between 
the components of X can be represented by means of a graph as follows: 




Xi -L Xj | Xy\{jj}. 


The graphical Gaussian models with edge and vertex sy mmetries, w hi ch we here call the 
colored graphical Gaussian model, have been introduced in IHoisgaard fc Lauritzenl [2008 ]. 
These models are defined as graphical Gaussian models with three different types of sym¬ 
metry constraints: equality of specified entries of the inverse of the covariance matrix K, 
equality of specified entries of the correlation matrix or equality of specified entries of K gen¬ 
erated by a subgroup of the automorphism group of G. These models are called respectively 
RCON, RCOR and RCOP models. In this paper, we consider only RCON models which 
form a natural exponential family with the precision matrix K as the canonical parameter. 
The model can be represented by colored graphs, where edges or vertices have the same 
coloring if the cor r esponding elem ents of the precision matrix are equal. 


Hpjsgaard & Lauritzenl [20081 ] proposed an algorithm to compute the maximum like¬ 


lihood estimates of K. However, to the best of our knowledge, there is no work for 
Bayesian est i mates o f K. Since the RCON model is a natural exponential family, we use 
the 


Diaconis fe Ylvisaker 


19791 ] (henceforth abbreviated DY) conjugate prior for K. This 


yields a distribution similar to the DY conjugate prior for graphical Gaussian models but 
with the symmetry constraints mentioned above. We will therefore call this conjugate prior 
the colored G-Wishart. 

Our sampling scheme is an adaptation of the independent Met ropolis - Hasti ngs ( 
forth abbreviated MH) algorithm for the G-Wishart proposed by 


Hence- 


Mitsakakis fe al. 


2011 ]. 


In the case of regular graphical Gaussian models (not colored), the DY conjugate prior for 
K is the so called G-Wish a rt and there are a number of sampling sch e mes for this distri¬ 


bution: see 


Lenkoski 


Piccioni 


200i 


1 


Mitsakakis fe al. 


2011 


Dobra fe al 


2011], Wang & Li 2012], 


20131 ], the more recent ones being generally more efficient than the preceding ones. 
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H owever, by t he v ery nature of these sampling schemes, only 


Mitsakakis fe al. 


Dobra fc al 


adapting 


20111] and 


2011 ] could be adapted to the colored G-Wishart. Moreover, we found that 


Dobra fe al. 


2011] leads to significant autocorrelation. The sampling scheme that 


we propose in Section 3 of this paper is therefore an adaptati on to the co lored G-Wishart of 


the sampling scheme for the G-Wishart given by 


Mitsakakis fe al. 


2011 ], 


In order to judge the accuracy of our sample, we need to either know the exact value of 
the expected value of the precision matrix or we need to proceed by simulations. We will do 
both. The RCON models for X = ( X v ,v £ V) that we consider in this paper are natural 
exponential families with density of the form 

f{X;K) oc exp{(IL, XX*) - ^\og\K\} 

where \K\ is the determinant of the precision matrix and [A, B) = tr AB denotes the inner 
product of the two symmetric matrices A and B. Let Q be the colored version of G to be 
defined in Section 2. The DY conjugate prior for K is then of the form 

P{K ; <5, D) = exp{—i«A- D) — (5 — 2) log \K\)} 

which is itself an exponential family in K and thus the expected value of K is given by 
the derivative of log Ig(S,D) with respect to the canonical parameter. The difficulty is, of 
course, to compute the normalizing constant Ig(5,D). 

Even for the uncolored G-Wishart, until recently, one did not know how to obtain the 
analytic expressi on of Iq(6,D) unless G was decomposable. Using an iterative method and 


special functions, 


Uhler et al 


2014] seem to have solved this very difficult problem but their 


results do not extend to our colored G-Wishart with symmetry constraints which, as we 
shall see, add another level of difficulty. We do give, in Section 4, the explicit analytic 
expression of the normalizing constant of the coloured G-Wishart for some special graphs Q: 
general trees, star graphs, a complete graph on 3 vertices and a simple decomposable model 
on 4 vertices with various symmetry constraints. For these particular coloured graphs, the 
analytic expression of Ig(6,D ) allows us to verify the accuracy of the estimate of the prior 
mean obtained with the sampling scheme developed in Section 3 by comparing it to the true 
prior mean obtained by differentiation of Ig{5, D ). The numerical results are given in Section 
5. In Section 6, we compute the posterior mean for two relatively high-dimensional examples 
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where the underlying graphs are cycles of length 20 and 30. Since we cannot compute the 
analytic expression of the normalizing constant in those cases, we compare the posterior 
mean estimate obtained through our sampling method of Section 3 with the true value of 
K used for simulating the initial Gaussian data. Our numerical results show good accuracy 
and small relative errors that decrease with sample size, as expected. 


2 Preliminaries 


We will now recall some definitions and concepts that we will need in the sequel. Let 
G = (V, E) be an undirected graph as defined in the introduction. Let Pq be the cone of 


p x y posi t ive c 


(see 


Lauritzen 


efinite matrices X with entry Xij = 0 whenever (i,j) 0 P- It is well-known 


19961 ]) that the graphical Gaussian model Markov with respect to G is the 


set of Gaussian 1V(0, £) distributions 

A/- g = {JV(0,£) | K = S" 1 eP G }. 


(1) 


The Diaconis-Yl visaker co njugate prior for the parameter K is the so-called G-Wishart 
distribution (see 


Roverato 


20021 ] 1 defined on P G and with density 


p(K\5,D) = —l—\K\^- 2 y 2 exp{-^(K,D)}, 

where 5 > 0 and D , a symmetric positive definite px p matrix, are the hyper parameters of 
the prior distribution on K and J G (<5, D) is the normalizing constant, namely, 


Ig(S,D)= I \K\( s ~ 2 y 2 exp{-l(K,D)}dK. 

Jp G 2 

Let us now define the RCON model. Let V = {14,..., 14} form a partition of V = {1,... ,p} 
and let £ = {Pi,..., E{\ form a partition of the edge set E. If all the vertices belonging to 
an element V) of V have the same colour, we say that V = {14,..., 14} is a colouring of V. 
Similarly if all the edges belonging to an element Ej of £ have the same colour, we say that 
£ is a colouring of the edges of G and that (V,£) is a coloured graph. 

Consider model {1]). If, for K € Pq, we impose the further restrictions that if 
(Ci) : m is a vertex class in V, then for all i € m, Ka are equal, 

(C 2 ) : s is an edge class in £, then for all (i, j) € s, the entries Kij of the precision matrix are 
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equal, then model (H|) becomes a coloured graphical Gaussian model called the RCON(V,£) 
model. 

For the computation of the analytic expression of Ig(5,D), we will need two special 
functions, the Bessel function of the third kind and the hypergeometric function p F q . The 
Bessel function of the third kind is defined as 

/*00 1 o 

K\(z) = / u 2X ~ 1 e~^ +u ] du. 

Jo 

For some special values of A, the Bessel function can be given explicitly 

K m (z) = K 3/2 (z ) = +z -3/V', 

K„ n {z) = ^0-' n + 32- S/2 + 

We will also use the classical formula 




K\ (y/pq) = I 

Jo 


, 2A - 1 o _ 5(^'+9“ 2 ) 


’du. 


The hyper geometric function p F q is defined by the power series: 

J? ( u u \ _ • • • i a p)k Z k 

p .F q (a u ...,apM,...,b q ,z)-'£ {h)k {bp ) k kl 


where 

1 if n = 0 

a(a + 1) ■ • • (a + n — 1) if n > 0 . 

The derivative of the hypergeometric function p F q (a\,... , a p , b\,. .., b q \ z) is given by 


(a)fc = 


ipFq ( a i > ■ • ■ > a p'i b) ■ ■ ■ > ")] — ^ ^ {pFq( a i +1, • • •, a p +1; b\ +1,..., b q +1; z)) . (2) 

3 The coloured G-Wishart distribution: a sampling method 

3.1 The coloured G-Wishart 

For G an undirected graph, let Q = (V,<?) denote its coloured version as defined in 
Section 2, and let Pg denote the cone of p x p positive definite matrices in Pq which also 
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obey the symmetry constraints of Q, i.e. 


Pg = { K G Pg | (Ci) and (G 2 ) are satisfied} . 


We define the GG-Wishart, i.e. the colored G-Wishart, to be the DY-conjugate prior for 
the parameter I\ of the RCON(V,£) model. Its density is 

p(K\8,D) = 7^^yl^l (<5 ' 2)/2ex P{-^ tr (KD)}l Pg (K) (3) 

where <5 > 0 and D, a symmetric p x p matrix, are hyper parameters and Ig(5,D) is the 
normalizing constant, namely, 


Ig(5,D)= [ exv{-\tr(KD)}dK. 

JPg 2 


(4) 


We will see that Ig(5,D) is finite only for D in the dual cone of Pg which has to be 
determined for each Q. We will derive the dual cones in the special cases that we consider 
in Section 4. 


In this section, following what has been done in 


Mitsakakis fc al. 


2011], we want to 


der i ve a M H algo rithm to sa mple from the GG-Wishart. But in order to do so, following 


Atav-Kavis & Massam 


200.5], we want to express the density of the GG-Wishart in terms 


of the Cholesky components of K scaled by D. To this end, we consider the Cholesky 
decomposition of ZW 1 and K written as 


D- 1 = Q t Q, K = 


with Q = (Qij)i<i<j< p and <h = {&ij)i<i<j<p upper triangular matrixes with real positive 
diagonal entries and we will use the variable 


T = 


-1 


We adapt the MH algorithm given in 


Mitsakakis fc al 


201 1] for the G-Wishart to the GG- 


Wishart, using the variable T rather than K. To do so, we first define 


v u (G) = min{(i, j) :i<j , (i,j) € u <E VUf} 
where the minimum is defined according to the lexicographical order and 

v(G)= U v u (G). 

uGVUS 
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We will write K V ^ G ' ) = (K%j | ( i,j ) € v(G)) for the free elements of K. The zero and colouring 
constraints on the elements of K determine the free entries : (i,j) € v(G)} and 

^,«(G) _ . (T j'j £ n(G)} of the matrices <f> and T respectively. Each non-free element <&ij 

and Tjj with (i,j) ^ v(G) is a function of the free elements and that precede 

it in the lexicographical order. The following two propositions give the expression of the 
non-free entries in function of the free ones and the free entries of K. The first part of each 
proposition can be found in 


Roverato 


2002 ], 


Proposition 3.1 Let K = <h T< I > be an element of Pg. Then the entries are such that 


&ij — 


2—1 

-K-ij ^2 ^ki^kj 
k =1 


k = o, 


= - 


<I>, 


2—1 

I] $ki$kj 
k =l 


®ij — 


in 1 2 1 

^iuju^i-uiu T X/ ^kiu^kju ~ S ^ki^kj 
k =1 k =1 


for (■ i,j) € v(G) 
forKi k = 0, k = 2,...,p 

(5) 

forKij = 0,... ,p, t / 1 
forKij^ 0, (i,j) € u € VU5, (i,j)£v(G) 


iu 1 


2 — 1 


= i<i> 2 . + 

w** ' / y ^ki u Z_j^kt I : 


( 6 ) 

/ori = 

(7) 

where ( i u ,j u ) = min{(i, j) : i < j and ( i,j) € u € V U £1} in t/ie lexicographical order. 


k =1 


fc=i 


Proof. The first three equations can be found in R overato [20021 ]. We will only prove ([6]) 
since © will follow immediately from it. For all (i,j) € u € Vu£ and (i, j) / ( i U i3u ) € u, by 

2-u 2 

©, we have that K iuju = ^ and in general K tj = ® k i$kj- Since K tj = K iuju , 

k= 1 k =1 

it follows that 

i u — l i —i 


^iui^iuiu + Y ®ki u ®k ju = $n$ij + Y 


kj- 


fc=l 


k =1 


Equations (|6|) and © follow then immediately. 
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Proposition 3.2 For K = Q T (fif T ^t)Q € Pg with 'I' and Q as defined above, the entries 
*1 \j of 'b are as follows: 

S — 1 

( ) CD 

( 8 ) 


Vrs = E^OTf 1 + ^ for {r,s) € v(G), r s , 

^CSS 


j=r 


$ 


'Ih, —- for (r,s) € v(G), r = s , 


s —1 


^rs = E 


j=r 
s — 1 


Q__ 

rJ Qss 


JS 


r-, + E -1 

E<- W - i l ' l ’-’-Y.' l, ',Tr ± ! f°rK„=0, rjtl, 


2—1 






*1. = E("^n ) forKu = 0, 

j =1 

i u —1 r—1 

+ X/ ^fciu^fcju X] ^fcr^fcs s —1 


'I'r-.s = 


fc=l 


k =1 


E^?r forK rs + {\ 


*&rrQ s 


j=r 


Q 


(9) 


(r, s) € u € V U £, (r, s) ^ u(G) , 


i—i 


l®L + eV- e*li» 


I 1 .. = 


k =1 


fc=l 


g s 


fors = 1,... ,p . 


( 10 ) 


Proof. We will prove Q and therefore (1101) . Since <b = for r fi s, we have 

S —1 

^rs — 'J PsQss X . ^ r jQj 
On the other hand, by ([8]) . we have 


cjs- 


j=r 


<i? r o — 


i u 1 v 1 

^iuju ^ + Z ^kiu^kju - Z ^fcr^fcs 

/c=l k =1 

^rr 


It then follows that 


s —1 


'i trsQss P X y ^rjQjs 


i u —1 r—1 

~b Z Z ^kr^ks 

k =1 fc=l 

$rr 


j=r 


which implies Q and (fT0|) . ■ 

Next we compute the Jacobian of the change of variable from K to 'ip v ( G " > in two steps. 
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Lemma 3.1 Let K be in Pg. Let vf be the number j €{*,... ,p} such that ( i,j ) 0 v(G). 
Then the Jacobian of the change of variable K v ^ —>• 3>^ G ) as defined above is 

det (J(K V ^ -> $^ G) )) = 2 |V| Yl 

i=1 

where |V| is t/ie number of vertex color class of Q. 


Proof. Order the elements of both matrices K and <f> according to the lexicographic order. 
For (i,j) € v(G), differentiating ((SJ) yields 


dK a _ ^ 
d$u * 


()K n 

d$ks 


= 0 for (fe, s) > (i, i) 


5ET 

—^ = 0 for (k,s) > (i,j), i / j. 
O^ks 


( 11 ) 

( 12 ) 


Therefore, the Jacobian is an upper-triangular matrix and its determinant is the product 
of the diagonal elements. The lemma then follows immediately from the fact that for i € 
{1 ,p} given, the cardinality of the set {(i, j) € v(G), ( i,j) > ( i, i )} is p — i + 1 — vf. ■ 


Lemma 3.2 Let K be in Pg. Let df =| {j : j < i, ( j,i ) ^ v(G)} |. The Jacobian of the 
change of variable (& V ( G ^ —>• vJ/HG) where and T are as defined above is 

det-> T^ g ))) = X\Q\f d? ■ 

1=1 

Proof. Order the elements of both matrices and T according to the lexicographic order. 
For (r, s) € v(G), differentiating ([5]), we obtain 


5$ 

d'F 


rs 


ss 


Q s 


d$rs 

d^ij 


0 for (i,j) > (r, s). 


(13) 


The Jacobian is thus an upper-triangular matrix and its determinant is the product of the 
diagonal elements. The lemma follows from the definition of df . m 


Theorem 3.1 Let Q = (V, £) be an arbitrary p-dimensional colored graph. Then the density 
of the CG-Wishart distribution expressed in terms of '& v ( G l is 


p(^ v ^\6,D) = 


2l v l 


V P 

p-vf-df+5-1 TT p- 


UT 


n*i 


—i—Vp+S—l 2 


i=l j=i 


(14) 


W,D)^ 

Proof. The expression of p(& v ( G fS, D) above follows immediately from the fact that \K\ = 
v 

n Tfj, that ( K,D) = fff-x anc l from the expressions of the Jacobians given in 

1=1 

Lemmas o and [3T2l ■ 
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3.2 The sampling algorithm 

We now briefly describe the MH algorithm we use to sample from the density (|14[) . We 
first note that if we make the further change of variables 

€ v(G),i 7 ^ j) •—> (t u = € v{G),^ ij ,(i,j) G v(G),i ± j) , 


we obtain 


p(tu,{i,i) e v{G),^ij,(i,j) e v{G),i ^ j\5,D) oc 


i 1 v o. -i V V 'll 2 . 

-1 2 ^ 2 ^ % J 


i= 1 j=i+l 


2—1 


5 - _I W 

- 1 2 1/11 O 

e i=1 has the form ofay . r r distribution. 

/v p— l—vY -\-0 


and we observe that f- 2 

We denote by \IfM and 'kt s+1 ] the current state of the chain and the next state of the 
chain, respectively. We denote 'k' the candidate for vp[■ s + 1 l. We a l so use the notation 


^v(cy = (*ij,(i,j) e«(G)j 

where v{G) c is the complement of v(G) in V x V. For (i. j) € v(G), an element 'I'|j is 
updated by sampling a value TF from a normal distribution with zero mean and standard 
deviation equal to one. For ( i,i ) G w(G), a element d'jy is updated by sampling a value 
('S’iiY from a chi-square distribution with p — i — vf + 8 degrees of freedom. The non-free 
elements of \E ,/ are uniquely defined by the functions in Proposition 13.11 and Proposition 13.21 
The Markov chain moves to ’F / with probability 


min{ 


h[(^')v(Gy] 

M(* M kG)c] 


1 }, 


where 


H^v(Gy) — n 


T p-i-v?+S—1 , 1 

1 exp(-- Y, 




2 z_/ ~ y > 

(■i,j)ev(G) c 


Finally, we can obtain iiM = Q T (\kM) T \E , M(5. Since h( x & v (G) c ) is uniforml y bounded b y L 
t he ch ain is uniformly ergodic (the strongest convergence rate in use, see lMengersen & Tweedie 
199f|). 


We now have a method to sample values of K from the CG-distribution, whether it is 
as a prior or a posterior distribution and thus obtain an estimate of the posterior mean of 
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K. In our MH algorithm, the candidates are drawn independently of the current samples 
through the proposal density. Thus, the algorithm gives an independence MH chain. Our 
simulation results in Section 5 will show that the chain has good mixing, low autocorrelation 
and high proximity to the true GG-distribution. 

The sample mean will converge to the expected value of K. In order to verify the accuracy 
of our sampling algorithm, we therefore would like to have the exact value of the expected 
value of K under the GG-Wishart. This is done in the next section for some special coloured 
graphs. 


4 The exact expected value of K in some special cases 


4.1 The mean of the CG-Wishart 


For a given Q, the GG-Wishart as defined in (|3|) and (JT]) clearly form a natural exponential 
family of the type 


f(K]9)dK = exp{(iG 9) — k(6)}fi(dK) 


with generating measure fi(dK) = \Ky- 5 ~ 2 ^/ 2 \p Q {K), 9 = — and cumulant generating 
function k(S,D ) = log Ig(5,D). To verify the accuracy of the sampling method given in 
Section 3, we will compare the expected value of K under the GG-Wishart and the sample 
mean obtained from a number of iterations of our MH algorithm. From the theory of natural 
exponential family, we know that the mean of the GG-Wishart is 


E(K) 


d k(5, D) 
d{~\D) 


n d k(8, D) 
dD 


We therefore need to determine for which values of 5 and D the quantity Ig(5,D ) is finite 
and then compute the analytic expression of Ig(5,D). We need also to differentiate this 
expression. 

We cannot do this in general but we will now consider several particular coloured graphs 
for which we can compute Ig(5,D). For the corresponding RCON models, we will see that 
when 5 > 0 (except in the case of the star graph with all leaves in the same colour class 
where we must have 5 > 1), the normalizing constant Ig(S,D ) is finite when D belongs to 
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the dual Pg of Pg. For any open convex cone C in R n , the dual of C is defined as 

C* = {yeR n \(x,y)>0, Vx € C \ { 0 }} 
where C denotes the closure of C. 

In the remainder of this section, for each RCON model, we determine Pg and the value 
of Ig(5, D ). This will allow us, in Section 5, to verify the accuracy of our sampling method. 

All proofs for Section 4 are given in Appendix 1 in the Supplementary hie. 

4.2 Trees with vertices of different colours and edges of the same colour 

Let V = {1,... ,p}. Let T = (V, E) be a tree with vertices of different colours and edges 
of the same colour. An example of such Q is given in Figure [Qa). Let a = (a*, i = 1,... ,pY 
where a* > 0 and b € R. Let S be the space of symmetric p x p matrices. We define the 
mapping 

m : (a, b ) € RP +1 m(a, b) € S (15) 

with m(a,b ) satisfying the conditions 

[m(a,b)]u = a,i, [m(a,b)\ij = b= [m{a 1 b)\j i for (i,j) € E , [m(a,b)]ij = 0 for (i,j) 0 E. 

Let M{Q) be the linear space of matrices m(a, b) for (a, b ) € R p+l . Let P be the cone of 
p x p symmetric positive definite matrices. Then 

Pg = M(G) n P . (16) 


Proposition 4.1 Let T be a tree as described above. The dual cone Pg is 

Pg = {m(a! , b') € M(Q) \ a' = (a'^i = 1 ,. .. ,p), b' € R, \b'\ < — j- ^ y^'a' } • (17) 

We are now in a position to give the analytic expression of Ig(5,D). 


Theorem 4.1 For Q = T as described above, 5 > 0 and D = m(a ', b') € Pg, the normalizing 
constant Ig(5,D) is finite and equal to 


/ e («5,zi) = 2l + p- 1 r(M[](a'; 


di—2 


5 

4 POO 


\i= 1 


I<s(\b\J^) \b\fe-^ W db (18) 

-°° ' (i,j)eE 
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where di denotes the number of neighbours of the vertex i in the tree (V., E ). 

For 5 = 1, we have 

/e(l,£>) = ( 2 vr)l Y ( a 'i a 'j) ~ (P ~ ^T 1 + i Y ( a i a i) + (P“ 1 ) b 'r i )- 

*=i (i,j)eE (i,j)eE 


For 5 = 3, Let er fc the kth elementary function of the variables yja^a'-, ( i,j) € E. We have 

Iq (3; D) = 2f” 1 7r§n(a')-f|>,r(fc + l)([ E (a'a')5-(p-l)6']- (fc+1) 

i =1 k =o (■ i,j)eE 

-[ £ (o'o;)i + (p-l)6']- ( ‘ +1) )- 

(i,j)£E 


4.3 The star graph with its n leaves in one colour class 


An example of star graph with its n leaves in one color class and different colors for the 
edges and the central node is given in Figure [jjb). For a € R,c € R,b = (fei, ..., b n ) € R n , 
let L{Q) the linear space of matrices of the form 


l(a, b , c) = 


a b\ b 2 

b\ c 0 

b 2 0 c 


u n 

o 

o 


0 0 


It is easy to see that the determinant of l(a, b , c) is 

\l(a,b,c)\=c n (a-^pj (19) 

and therefore, Pg is the open cone 

llftll 2 

Pg = {l(a, 6 , c) € L(G) : c > 0, a - — > 0.} 

c 

The dual cone P*{G) and the normalizing constant Ig(5,D) are given below. 

Proposition 4.2 For a star graph with all n leaves in one colour class, the dual of Pg is 

P* = {l( a ',b',c') € L(G) | ||fe'|| 2 < na'd}. (20) 
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Theorem 4.2 For Q a star graph with all n leaves in the same colour class, 5 > 1 and 
D = l(a' ,b' ,d) € Pg, the normalizing constant of the CG-Wishart is 

I S (6,D) = X rfCi-xx—x> x (nn , c , H|t ,| |2) „_ 1)t+1 X w- 1) = + 1)1$. 


4.4 The star graph with all vertices in one colour class 


An example of star graph with all vertices in one color class and different colors for the 
edges is given in Figure (He). This case is a special case of the preceding one and therefore, 
we have immediately that 

Pg = {l(a,b,a) € L{Q) \ a > 0 , a 2 — 116|| 2 > 0 }. 


Since this is a well-known cone, called the Lorentz cone, we know also that it is self dual 
and therefore 

Pg = {l( a ',b',a') € L(Q) \ a' > 0 , {a'f - \\b'\\ 2 > 0 }. 

It remains to compute Ig(5,D). 


Theorem 4.3 For Q the star graph with n leaves and all vertices in the same colour class, 
5 > 0 and D = l(a',b',af) € Pg, the normalizing constant of the CG-Wishart is 


(n+l)5 

> O -t 


Ig(S,D) = 


where u = (and 2?(|, |) is the Beta function with argument (|, § 


C„r((n + !)§)„ .8 n 


(re+l)<5 


(re+l)5 


(n-fl) 2 (o') 2 


'2 2 ' 


-B( —, o) 2Pi (( n + 1)7, {n + 1)7 + 


1 n + 6 
2 ’ 2 ’ 


4.5 A complete graph on three vertices with two edges in the same colour 
class 

This graph is represented in Figure QJd). In this case, the cone Pg is the set of positive 
definite matrices K = (&r/)i<ij <3 with &13 = fc 23 . 

The dual cone P*(Q) and the normalizing constant Ig(5,D ) are given below. 

Proposition 4.3 For the graph in Figure\l\(d), the dual of Pg is 

Pg = {P = {dij)l<i,j<3 € S | di3 = ^23) 

da > 0, i = 1 , 2 ,3, d \ 2 < ^ 11 ^ 22 , 4 df 3 < (dii + e ?22 + 2di 2 )d33}- 
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Figure 1: Black vertices and edges are all of different colours, (a) The colored tree, (b) 
The coloured star with the centre vertex of a different colour, (c) The coloured star with 
all vertices of the same colour, (d) The triangle with two edges of the same colour, (e) The 
decomposable graph with three different colours for the edges. 

Theorem 4.4 For Q as in Figure\l\(d), 5 > 0 and D G Pg, the normalizing constant of the 
CG-Wishart is 

Ig(S,D) = 2 2 vrr(-)(r(—-—(dn + c ?22 + 2 di 2 ) 2 [^ 33(^11 + d 22 + 2 di 2 ) — 4df 3 ] ~ 

x{dnd22 — d 12 ) 2 ■ 

4.6 A decomposable graph with three vertex classes and three edge classes 

This graph is represented in Figure [He). Then the cone Pg is the set of matrices of the 
form 

^ ku ki2 k 13 h 4 ^ 

ki2 k 22 k\ 3 ki4 

K = 

ki3 k 13 k 33 0 

\ ku ku 0 k 33 ) 

Proposition 4.4 For Q as in Figure^ e), the dual cone is the set of matrices 

Pg = {D = (djj ) \<i.j <4 € S | d 23 = di 3 , d 2 4 = di4, ^44 = d 3 3 , d\\ > 0 , 

dnd 2 2 — d\ 2 > 0, dn + 2di2 + d 2 2 > 0, d 3 3(dn + 2di2 + d 2 2) — 2(di 3 + ^ 14 ) > 0} • 
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Theorem 4.5 For Q as in Figure^ e), 5 > 0 and D E Pg, the normalizing constant of the 
CG-Wishart is 

Ig(6,D) = 2 s+ \lr( d -)r(^)r(5)(d 11 +d 22 + 2d 12 ) 5 - 1 (dnd 2 2-d 2 12 )-^ 1 
x [<^33(^11 + d 22 + 2^12) — 2 ( 4^3 + d 2 4)] <5 . 


5 Numerical experiments when we know the exact mean 


In order to illustrate the performance of our MH algorithm, we conduct a numerical 
experiment for each of the colored graph (a) - (e) shown in Figure [TJ In each case, for 
a given D and S, we first derive log Iq(8,D), then the prior mean E(K ) under the CG- 
Wishart by differentiating logic(5, D) with respect to — y. We then sample from the CG- 
Wishart. We run the chain for 5000 iterations and discard the first 1000 samples as burn 

/n a y-\i=5000 

in. Our estimate K for K is the average K = ^’^ooo — 1 the remaining 4000 iterations 
Ki,i = 1001,... , 5000. For arbitrary K and K' we define the normalized mean square error 
( nmse ) between K and K' to be 


nmse(K, K') = 


\K-K’\\l 

\\ k '\\ 2 2 


where | |ii| || is the sum of the squares of the entries of K. We repeat the previous experiment 
100 times, obtain K J , j = 1,..., 100 and compute 


1 100 

Wmse(E(K),K) = —J2nmse{E(K),K j ) 

3 = 1 

where E(K) is obtained by differentiation of log/g(<5, D) with respect to — y at our given 
D and 5. 

For each graph in Figure [T] for an arbitrary j £ {1,... , 100}, we give the trace plot of 
log | K'l |, i = 1000,... , 5000. The traceplot shows that the chain seems to be mixing well. 
We also provide the autocorrelation plot with time-lag h for log \Kf\, i = 1000,... , 5000 in 
function of h where, for an arbitrary given j , we define the autocorrelation coefficient for 
Yi = log | Kf |, i = 1000,..., 5000 to be 


Rh 


EZTooo (Yj-YKYj+b-Y) 

V 5000 (Y■ — Y ) 2 

iCi=iooo G* 1 ) 
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Q 

<5 

log Ig(S,D) 

nmse{E(K), I \) 

Fig. DJa) 

1 

7 

~\ E log a'+ log 
2—1 

i i 

Efa'a'p—66' E(a'a')2+6 V 

-i= 1 i= 1 


0.0069 

Fig. ETb) 

3 

^loga' — 91og(8a'c' — \\b' 2 ) 

0.0187 

Fig. HXc) 

3 

15 log a' + log 2 -F 1 ( 2 5 ,8; 6; 

0.0064 

Fig. Did) 

3 

| log d - 2 log(d 33 d - 4df 3 ) - 2 log(dnd 2 2 - d\ 2 ) 

0.0005 

Fig. QJe) 

3 

2 log d — 3 log(d 33 d - 2 d\ 3 - 2 d\ A ) - 2 log(dnd 2 2 - d\ 2 ) 

0.0009 


Table 1: For the graphs of Fig. Q] and 5 given: analytic expression of log Ig(5, D) where 
d = dn + ^22 + 2di2 and value of nmse(E(K), K) averaged over 100 experiments. 

The autocorrelation plots indicate that the samples have a low autocorrelation. The nu¬ 
merical values of the matrices D, E(K) and K as well as the traceplot and autocorrelation 
plot of log (| K |) for all five graphs in Figure [T] are given in Appendix 2 in the Supplementary 
file. An overview of calculations and results are given in Table Q] which, for all different 
five colored graphs in Figure [T] shows the parameter 5 we chose for the prior distribution, 
log Ig(5,D) and the normalized mean square errors. In order to obtain the mean E(K) of 
the C'G'-Wishart for the graph in Figure (He), we use formula (|2|) to get the derivative of 
the hyper geometric function p F q (ai ,..., a p ; bi ,..., b q ; z). We see that the normalized mean 
square error is of the order of 10~ 3 or less except for the star graph with all leaves of the 
same colour in Figure (HT). 
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(a) 


(b) 


(c) 


Figure 2: Cycles of length 6 with the three different patterns of colouring that we use for 
the cycles of length p = 20 and p = 30. Black vertices or edges indicate different arbitrary 
colours. 


6 The posterior mean from simulated data: p = 20, p = 30 

In this section, in order to assess the accuracy of our sampling method for larger graphs, 
we generate data from a N(0, K -1 ) distribution with K given in Pg. We take the CG- 
Wishart with 5 = 3 and D = I as the prior distribution of K. Clearly the posterior 
distribution will be C'G'-Wishart with parameters 5 + n and / + nS where S is the sample 
covariance matrix. We will use this posterior and our sampling method of Section 3 to 
compute the posterior mean E(K\S) as an estimate of K. 

We run our experiment with six different coloured graphs. For three of them, the skeleton 
is a cycle of order p = 20 and for the other three, the skeleton is a cycle of order p = 30. 
For each cycle of order p , we give three different patterns of colouring which, for the sake 
of saving space, are illustrated in Figure [2] for p = 6. The values for the entries of K for all 
three types of graphs are as follows: 

Ku = 0.1, i = 1,3,... ,2p- 1, Ku = 0.03. i = 2,4,... ,p, 

I<i,i +1 = K i+1}i = 0.01, i = 1, 2,... ,p - 1, Ki p = K p i = 0.01. 

Though, for convenience, we chose, for all the models, the same values for the entries K{j,i ^ 
j to be all equal to .01 and for K tl to have two different values .1 and .03, in our computations, 
we used, of course, in each case, the model represented by each of the respective graphs. 
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For each graph, we generated 100 datasets from the N(0, AT -1 ) distribution. The posterior 
mean estimates are based on 5000 iterations after the first 1000 burn-in iterations. We 
denote K = (Kij)i<i,j<p the posterior mean estimate. 

Table [2] shows nmse(K , K ) for the three colored models on the simulated examples 
when p = 20 and p = 30, averaged over 100 simulations. Standard errors are indicated 
in parentheses. Computations were performed on a 2 core 4 threads with i5-4200U, 2.3 
GHZ chips and 8GB of RAM, running on Windows 8. We also give in Table [2] the average 
computing time per simulation in minutes. 

Table 2: nmse{K, K) for the three colored models when p = 20 and p = 30. 


g 

p = 

nmse(K, K ) 

20 

Time/sim 

p = 

nmse(K, K ) 

30 

Time/sim 

Fig-2 (a) 

0.005 (0.003) 

19.425 

0.040 (0.021) 

86.423 

Fig.2 (b) 

0.011 (0.003) 

18.739 

0.033 (0.011) 

82.876 

Fig.2 (c) 

0.039 (0.021) 

16.410 

0.080 (0.033) 

82.563 


In Table [31 for the graph of Fig. [2] (a) with p = 20 and p = 30, we give the values of 
the entries of K together with their batch standard errors. For the other models, average 

Table 3: The average estimates and batch standard errors for K in Fig.2 (a). 


P 

K u 

K V2 

K lp 

P‘2'2 

20 

0.1040 (0.0005) 

0.0103 (0.0002) 

0.0104 (0.0002) 

0.0313 (0.0001) 

30 

0.1223 (0.0009) 

0.0121 (0.0004) 

0.0125 (0.0004) 

0.0361 (0.0003) 


values of the entries together with batch standard errors are given in Appendix 3 in the 
Supplementary file. 

Remark 1. At this point, we ought to make an important remark. In Section 5, we 
proved that the CG-Wishart was proper for D € Pg. When we compute the posterior 
mean in this section or more generally for any colored graph, even if D > 0 belongs to Pg, 
the hyperparameter D + nS does not usually belong to Pg of course and yet the integral 
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Ig(5 + n, D + nS ) converges. This is due to the fact that we can write S as 

S = Si + S 2 

where Si is the projection of S on the subspace of p x p matrices with fixed zeros according 
to Q and equal entries for edges and vertices in the same colour class and S 2 belongs to its 
orthogonal complement. Since K £ Pg, we have 

0 < (K, S) = (K, S{) 

and, since the inequality above is true for any K £ Pg, it follows that Si belongs to Pg. It 
follows also that Ig(S + n, D + nS ) is finite. 

Remark 2. For the computation of the posterior mean following our sampling scheme of 
Section 3, we may wonder whether we should take Q to be such that Q t Q = (D + nS)^ 1 or 
Q t Q = (D + nSi) -1 . We take Q to be such that Q f Q = (D+nS)^ 1 to use all the information 
given by the data. 


Appendix 1 


Proofs of Section 4 

Proof of Proposition 14.11 

Let M be the set of p x p matrices. Let 

7g = {x € M I Xij = 0, for i < j, Xij = Sij ^ 0, for i > j, (■ i,j) £ E, Xu = ti > 0, i = 1,. .. ,p} 


be the set of upper triangular matrices with positive diagonal elements and nonzero entries 
. i > j only for (i,j) £ E. The vector s = ( Sij , (i, j) £ E) belongs to R p_1 since a tree 
wit h y vertices ha s y — 1 ed ges and t = (U, i = 1 belongs to R p . It is well-known 

(see 


Paulsen et al 


19891 ] and 


Roverato 


2000]) that we can find a perfect elimination scheme 


enumeration of the vertices of T such that, with this enumeration, K £ Pg can be written 
as K = X(t, s) T X(t , s ) with X(t , s) £ 7 g- Then for K = K(a, b ) as in dTH we have 

a j — tj T ^ ' Sjj, b — 

ieEi 
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where (t,s) is the Cholesky parametrization of K 6 Pg. We can also parametrize K £ Pg 
with (t, b ) € (0, +oo) p x R using 



i£Ei i 


( 21 ) 


In this proof and the following one, we assume that the numbering of the vertices of T follows 


a perfect elimination scheme ordering. We then say that the last vertex p in that ordering 
is the root of the tree and we will write 


Ej = {i,i<j \ ( i,j) € E}. 

For convenience, we denote by C the right-hand side of equation (1171) . 

We show first that Pg C C. Let D = m(a' , b') € Pg. Using d2T]) . we have 

(K, D) = a i a\ + • • • + a p a p + 2 (p — 1 )bb' = t p a p + Ab 2 + 2 Bb + C > 0, (22) 


where 



(23) 


Now observe that for fixed i = 1,... ,p then either i = p and the set {j;i € Ej} is empty 
since p is the root of the tree, or the set { j : i € Ej} is reduced to one point, say jj. Therefore 
we have a 'j = a ii for i < p and zero for i = p. (For the graph in Figure [U (a), we have 

h = h = 35 = k = 7 and j 3 = j 4 = 6) and it follows that 


r >— 1 



(24) 


Let us prove that a' > 0 for all j = 1,... ,p. Take (ai,..., a p ) € [0, oo) p \{0,... , 0}. Then 
K(a,0) G Pg\{0} and (K,D) = a\a\+- ■ ■ + a p a! p > 0 implies that a '• > 0 for all j. Let us now 
prove that Ab 2 + 2 Bb + C > 0 for all b. If not, there exists bo such that Ab q + 2Bbo + C < 0. 
Since a' p > 0 taking t p very small and b = bo in (f22l) gives a contradiction. 

Let us prove that 





(25) 
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Since V 6 , Ab 2 + 2Bb + C > 0, we have B 2 < AC. Now consider the function 


(ti,..., t p - 1 ) AC 

and let us compute its minimum A*C* on (0, oo) p . This function AC is homogeneous of 
degree 0 and therefore if its minimum is reached at t* = (t*,, t*_i) it will also be reached 
on nt* for any k > 0. We have for i = 1,. .. ,p — 1 

AC = 2t i a' i A - \a':C = 0 

L * 4-6 Jl 


dt 


and we therefore have 


ti = K 


,/ X 1/4 

'0i_ 

J 


v-1 


. a - = c- = ■£ = E 


i =1 


a'a'. 


Since B 2 < AC for all (t\,... ,tp- 1 ) € (0, oo) p - 1 , we can claim that B 2 < A*C* or equiva¬ 
lently (1251) . 

Let us prove that inequality (1251) is strict, that is B 2 = A*C* is impossible. Suppose 
that B 2 = A*C*, i.e. | 6 '| = A*/(p — 1) > 0. Then with ti = t* we get Ab 2 + 2 Bb + C = 
A*( 6 +sign 6 ') 2 . Taking b = —sign b' and U = t*,i = 1,... ,p— 1 yields Ab 2 +2Bb+C = 0. Now, 
letting also t p = 0 in (1221) . we see that the left hand side of (122]) is zero for an (a, b) € S\ {0} 
which is not zero, since b = ±1. But this cannot happen for D(a',b') € Pg. Therefore 
is strict and the proof of Pg C C is completed. 


Let us now show that C C Pg. For D(a',b r ) € C given, we want to show that (K,D) 
is positive for all K(a,b) € Pg \ {0}. We will do so first for K(a,b) € Pg and then for 
K(a, b) € Pg \ {Pg U {0}). For K(a, b) € Pg, tk > 0 and b € R. From (122]) . we have 

( 6 + l) + ^ AC ~ B ^ • 

We have checked above that AC — B 2 > 0. Moreover a' p > 0 since D(a', b') € C. It follows 
immediately that ( K , D) > 0. 

Let us now show ( K , D) > 0 for K{a, b) G Pg \ {Pg U {0}) that is for t\ ... t p = 0 and 
(ti,..., t p , b) 7 ^ 0. We need only show that then ( K , D) ^ 0. But 0 = aia/| + • • • + a p a' p + 2{p — 
1 )bb' = Yl P i =i tja'i implies that t t = 0 for all / = l.... ,p since {a 1 , b') € C implies a! i > 0. But 


(K, D) = t 2 a p + Ab 2 + 2Bb + C = t 2 a' p + A 
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since b = Usij, this implies 6 = 0 but this is impossible since we exclude the zero matrix for 
K. 


Proof of Theorem 14.11 

In Ig(6, D ) we make the change of variable (12111 . Switching to these Cholesky coordinates 
leads to the Jacobian dadb = 2 p t\... t p dbdt. As seen before the new domain of integration 
is the product 

{( 6 , t)]tk > 0,b € R} = (0, oo) p x R. 

With the notation A, B, C of (1231) . we have 

(A (a, 6 ), D{a! , 6 / )) = 2 (jp — 1)66 / + + • • • + CLpCi'p = t p dp T Ab^ + 2Bb + C. 

Using (1241) for the expression of A, we obtain 


t 2 / P~ 1 t i a i 1,2 a jj 

Ig(8,D) = 2 p I (ti... t p ) s ~ 1 e~^ p ~ 1 ' >bb e ^ TT e 2 2t ; dt\... dt p db 

I (0, oo)PxR 


= 2 P 


'L 

l 


i =1 


v -1 


u p^p 

e — 


t S p 1 dt p f e (p 1)6b ' Yl fe /2 (|6|(a'a'.) 1/2 )(|6| Ja'./ai) 5/2 ) db 
J-oo j=i 

<*> 


with the notation 


MD) = j 

We now prove by induction that 


U=1 


oo pl 

-(p-l)bb'\ h \(p-l)5/2 




(27) 


2=1 


p— 1 / 

r n' 


K) 


* n £=n« 


d 7 ;-2 


^2 n 
2=1 2=1 


(28) 


Of course PSI) is correct for p = 2. Suppose that (|2S1) is true for any rooted tree with size 
p. Consider a rooted tree T* with vertices {0,1 ,,p} and root p and numbered, as usual, 
such that i -< j implies i < j. Denote T the induced tree with vertices {1,... ,p}. Finally 
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denote d* = (dg, ■ ■ ■, d*) and d = (di ,..., d p ) the number of neighbours in T* and T. Then 
dg = 1, d* 0 = 1 + dj 0 and d* = di if i / 0 and i ^ jo- This implies that 


1 


P~ 1 / 

t— r ^7 

xfl * 

P J i=0 


4 


(°p) 2 f n a i a 0 («p) 2 a, 


TT 4 

X 1I 


a,. (1) o' 


/ p p 

do ~X~\rJ\di—2 (?) 


P 1=1 1=1 1=0 

where (1) comes from the induction hypothesis and (2) from the link between d and d* . The 
induction hypothesis is extended to p + 1 and (PHI) is proved. 

We now prove that Js(D) defined by ([271) converges if D = m(a',U) € Pg where Pg is 
the convex cone defined in Proposition 3. We write Js{D) as the sum 


T“T- 2 =nw 


d*- 2 


Js(D)= [ ...db+ [ 

J —oo J 0 

When b —>• ±oo, |6| —>• +oo. From 


+oo 

db+ I ...db. 

—oo J 0 


(29) 


Watson 


19951 ) page 202, 7.23 (1) we have 


K x (s) 


7 r e 


2 s 1 / 2 ' 


We use this fact to analyse the convergence of Js(D). If D = m(a',b') € Pg, from the 
asymptotic formula above, we see that the integrands in both integral on the RHS of (1291) . 
when |6| goes to infinity, behave like \b\ c e~^ H where, since m(a! ,b') € Pg, 

p - - 

H = V 1 a i a j ~(p~ 1 )|& , |sign (bb') > 0 




and c = (p — 1)^-. Since the argument of (f27T) is continuous, both integrals converge at 
infinity. 

To study the convergence of these integrals when b —>• 0, we recall that 

r+oo x 

2K x {s) = / X +Pdx. 

Jo 

Making the change of variable u = sx in the expression of 2K\(s) we see that 


■'s-s-O 


s _A 2 A-i r ( A ). 


K x (s) 


Therefore, for both integrals in the RHS of (l29|) . the integrand is equivalent to 


(|6r 4 r(|)) P_1 |6|¥ e -<P-i)«»' = | 6 |f e 


— (p—l)bb' 
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and therefore both integrals converge at 0. The expression (|18[) of the normalizing constant 
is now proved. 


By (USD, I G (1,D) = 2^r(±)(a;)4(nf- 1 1 where 

p -1 


r 

Ji (D) = / 

J — < 




I^II KimMa'jf^db 


2=1 

P-1 


/ OO P 1 / 1 

e _ Cp _ 1 ) N ,.| 6 |E=l(-Q 

-°° i=1 v 

,7T,tl TT / , , ,_A Z -00 — (P-l)bb'-\b\ P E(aJaJ)4 

= ( ? ) 2 n^4) 4 / e i=i db - 

z i=i • y -°° 


We compute the integral 


/ 


00 -(p-l)66'-|6| „ 

e i=i db = 


L 


0 -(p-l)M>'+6 ^W.)^ /-°° -(p-l)6b'-b eV^'X 

e <=1 1 db+ e *=* 1 

Jo 

1 1 
-+ 


E («) - (p - l)*' E «<) + (p - !)&' 


2=1 


2=1 


Therefore 


, 7T, P -1 


P-1 


. P-1 


= {^Y 2 ~ n( a <°ji) 4 [(xxxx - (p- 1 ) &/ ) + (xxx-x+(p- 1 ) &/ ) 

«=1 i=l i=l 


P-1 


- 1 - 


P-l 

Since £ (a'a',) = E(ij) e sK a j)> this yields the expression of /g(l, ZD). 
i=l 


db 


25 





Similarly, from (JUJ), Ig( 3 ,D) = 2 p+ 2 r(|)(ap) 2 [I!=i ("^) *h(D) with 
/■oo , P -1 

Jsp) = / e -(p-D** | 6 | 3 Cp-D n^dftKa'a'J^)^ 

■ y -oo i=l 2 

/ oo rzr 1 

e -(P-i)»|i,|!(p-i) J][ / (|6|-iWaj,) -1 + |6r*(o'4)-J)e" l< ’ IW “ S < )3 ]<ii> 

*=1 V ^ 

/■OO P- 1 P- 1 1 

= / e- (p - 1, “>l* (p - 1, (|)' ?1 |(>r l<p - 1 ) IL^O" 5 nKK'IWO 4 + l)e‘ |i|K “i<P]<i6 

J -°° i =1 i=l 

= (^) 2 n^s) 4 / e i=i na+w^*) 2 )* 

t=l ' y_0 ° i=l 


T-r 1 ,//,-» /■“ -(P-I)W-IH E(«;«',)» 

(ji^IlW) j y e 


P-1 


1 '‘ (1 + |6|cri + |6| 2 fj 2 + ... + \b\ p l (j p -i)db 


i =1 


ir.ei?-!, , , ^4 /■“ -(p-i)W-IH E K«;,) J 

5 )’ nw4)‘E"‘ e 

Z ■ 1 1 J—OO 


\b\ k db , 


where the ct,; = <7j(^/a'a'-., i = 1,...,_p — 1) are the symmetric functions of yja^a'j., i = 
1,... ,p — 1. Since 


/ 
then 


— (p— l)bb'— 16| 5]) ( a i a i-) 2 1 

e *=1 * I6I"*- 1 * = r(m) 


p-i 


/g(3,D) = 2*" 1 tt* n (a i)“ 5 Z)^r(fc + 1) 

i= 1 k =0 

This yields the expression of Ig(3,D). 


p-i 


p- 1 \ _m /p-i 

^(a'a'.)5 - (p- 1)6' - I ^(a'a'J5 + (p - 1)6' 

i=l / \i=l 


P-1 


^(a'a'J 2 - (p - 1)6' - ( ^(a'a'Js + (p - 1)6' 


Proof of Proposition 14.21 

By definition Pg = {D = l(a',b',c') € M{Q) \ {K,D) >0, K G Pg\ {0}}. Let f3 denote the 
angle between b and b'. Then, since cos/3 > — 1 

(. K,D) = aa' + ncd + 2||6||||6'|| cos (3 > aa' + ncd — 2||6||||6 , ||. 

Therefore 2||6||||6'|| < aa' + ncd and since ac > 0, W < ^ aa + ™~ c ^ . By differentiation 

with respect to a and c, we see that yaa + ^ c ^ > Ana'c' and therefore (K, D) > 0 implies 
that H^H 2 < na'c'. 
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Proof of Theorem |4.21 Let 

us introduce the matrix 




r 

si 

S2 • • • 




0 

t 

0 ... 

0 

A(r, 

s,i) = 

0 

0 

t ... 

0 



0 

0 

0 ... 

t 


If (a, 6, c) G Pg the only triple (r, s, i) such that t > 0 and r > 0 and such that 
K(a, 6, c) = A(r, s, t)A T (r, s, t ) = 


r 2 + ||s|| 2 s'i 


is 


t 2 /n 


satisfies r = (a — ■l^-) 1 / 2 , t = -^/c, s = -^=. A new parameterization of Pg is therefore given 
by the change of variables ( a,b,c ) into ( r,s,t ) with a = r 2 + ||s|| 2 ,6 = is, c = i 2 , where 
(r, s, i) belongs to 


{(r, s,t);r > 0, s G P n ,i > 0} = (0, oo) x R n x (0, oo). 

With this parameterization, from (1191) . we have detPT = r 2 i 2n and dadbdc = 4 rt n+1 drdsdt. 
Then 


Ig(5,D) 



-\\s\\ 2 a'-2t(s,b') 


/ 0 J 0 




5- ds ) r s ~H^ n+1 e ~ ra r tc drdt 


7T \ n /2 r c 

a '' Jo 


e 2 


- t (5-l)n+l dt x 


poo 

/ r 5 - 1 

Jo 


n/2 


f «iih , n 2 


I —nvc | Q||Q || / c_-. \ n 

e 2 5^— v y ’ 2 dv x 

o Jo 


r 4-i 

/ t>2 

Jo 


e 2 dr 

~ a/ 

e 2 dv 


2 ^±2 7r n/2 a ,(f-l)(n-l) 


1 


77 

(na'd - ||6'||2)(^-i)f+i r(( ' 5 “ ^2 + 1)r( 2 )- 


Proof of Theorem 14.31 



Ig(5,D) = / I / a n 2 ( a — 1LLL) 2 exp— {{n + l)aa! + 2{b 1 b')}da\ db 

lR n \7||6|| a 2 y 

a ( ' Tl_1 ' ) “(a 2 — ||6|| 2 )^"~ exp —-{(n + 1 )aa! + 2(6, b')}da ) db 

hii 2 
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Let us make the change of variable 


(a, b ) € (||b||, +oo) x R n >-)■ (u, R , 0) € (0,1) x (0, +oo) x S 

R 

y/u' 


where b = R6 and S is the unit sphere in R n and a = -^=. We have dadb = — 2 Jjy 2 RC n R n 1 dadRdO 
where C„ is the surface area of S. Then 


'«<“■> - ^I.L 


R K ‘ 


5-2 


R S 2 (-1) 2 exp — { 


5-2 


(n + l)i?a / 


2y/u 


+ K<6», b')}u _3/2 du j R n dR 


c n r r+°° 

2 Js |/0 \j 0 

Cn / / + °° ( f 1 - R(n_1) ^^ 14_(n+1) ^ T ^- R,5_2 ( 1 - exp -{ — + -R(0, 6')}u _3/2 <iuj H n di? 


2 Js 


C n 

2 


/„ [r (/: 


„(n + l)|-l — (n + 1) ——j— — ,^ 5 — „ r > 

it v ' 2 U 4 2(1 — U) 2 exp— 


<5 — 2 (n + l)a / 


2y/u 


+ (0,6 )}du J dR 


dG 


C»r((n+ 1 ) j) 


2 y/u 


2 y/u 

■ + (0, v "h ^ 2 du 


Cn r f 1 (n I 11 g ~ 2 3 <5-2 / /*+oo £. _i (n + l)a , \ 

- / / t/, ^ 3 ^ (1 — it) 2 j / ^ exp —-R{- + (0,6 ) }di? ] du 

2 Js Jo \J 0 2^/57 / 

/■ /-i . i n 8-2 3 5-2 , 

■Is L u 


Cn r ((n + 1)^) (n + l)a/ _ f C ^ _f n xi^ <5 2 _3^ 5 2 (nA-1\ — / 2 , \ - 

2 2 ' >2 / s [/„ “ C + ) 4 2(1 - u) 2 “ ( + ), ( 1 + (^TTv^< 9i6> ) 

*»,«(<*') f f 1 u f-hi-u)!- 1 f^(-i) fc ( 2<e,i> } 

JS ^0 v (n + l)a/ 7 

(n+l)5 , 

■>-T-J- / 


2(6,6') yk k + 


where K n ^{a!) = 


2—T 


_ Cnr((n+1)^) rpj^ £ 

(n+l)S (n+l )5 • rncicioie 
(n+1) 2 (a') 2 


j: 


k =0 

oo 

= X n>5 (a')V( 


(n + l)a' 

2 N l 2fc(( n + l)|) 


tt~ -i (l-u)3 -i dit / (6,b') k d6 

Js 


t n» + i)«' 


u 2 —/ {0,bY k d6 

Js 


= K n , s {a') V( 


k=0 


(n + l)o' 


2t ((» + 1)1)^ r(i + f )r(f) 2t (i/ 2 ) t 


(26)! 

Q+l\ 


r(A: + ^) (n/2) fc 


We now use the fact that (a+fc = 2 2fc ( ^ I ( ++ j anc [ p( a _|_ — r(a)(a)fe. We also use 

\ z / k\ ~ / k 

the fact that 

(2fc)! = (135...(2fc - l))(246...2fc) = 2 k fc!2 fc i|...^i = 2 2fc fc!|(| + l)(i + 2)...(| + (fc - 1)) = 2 2fc fel(i) fc . 


Finally, since the integral is rotational symmetric, we take b' = ||i/||ei so that (0, (/) = ^i||6 / || 
and recalling that d6 is the distribution of tt5tt when Z ~ Af(0,1) so that 9\ = .= Zl = 

II II Y +...+-^n 

which is then such that Of ~ Beta(|, n ^-), for v = 0\, we have 


{0,b') 2k d0 = 


,/||2 k 




1 71—1 


2(1 — v) " 2 1 Mu = 


/1| 2fc 


( 1 / 2 )* 













































Writing B(a,(3) for the Beta function with argument (a,/3), we obtain 

S 


I^,D) = ((^ 7 ) (S)l(<" + l)4),((" + 1 )4 + 2) 


5 , h Wk „,e, 2 fc ( 1 / 2 )fc 
" (n/2) fc 


(^) 


= 


k 

,,,J n 5 J In (?) fc „ L /„2fc (1/2) fc 

,«(«)- B ( 2 ’ 2 ^S((n+l)a') 2 2k k\(^)t v n ^4 + 2 )*, ^ ^ (n/2) fe ' 

Let u = () • We note that since A = l(a',b',a') € Pg, then n < 1. After obvious 
simplifications in the expression above, we have 

A n 00 n,k ((n + 1)^1 f(n + 1)| + tj') 

w,d) = ~ h 


k =0 


/n+5\ 

V 2 7/ 


, /ln 7 n, ^ , ,>() , , d in + d \ 

= K nt s(a)B(-,-) 2 Pi(jn + l)-,(n + l)- + -,^—;oj 


5 l n + 6 


Proof of Proposition 14.31 

We write the Cholesky decomposition of K under the form K = AA l with 



( 

On 

012 

013 ^ 

A = 


0 

022 

023 


V 

0 

0 

033 / 


Expressing the kij in terms of the aij and imposing A 43 = k 2 z immediately shows that we 
must have 033 = < 223 . Then, let D = (c%)i<jj <3 with ^33 = ^23 since the dual of Pg must 
be in the same linear space as Pg. 

(A, D) = ( 0^3 + <232 + 033)^11 + ( a 22 T ®13)^22 + <^33^33 + 2(a22®12 + 033)^32 + 4033033(^33 

= af 3 (dn + d 2 2 + 2^32) + 4033033(^33 + a\ 2 du + 2 ai 2 a 22 di 2 + 033^33 + 0^22 + 033^33 , 


which we view as a quadratic form a l Ma with a 4 = 


M 


^ dll + d 2 2 + 2^32 2^33 

2(^33 ^33 

0 0 

0 0 

\ 0 0 


(oi 3 > a 33 ; a i 2 > 022, 033) and 

0 0 0 \ 

0 0 0 

d\\ di 2 0 

d\2 d 22 0 

0 0 0?3 3 J 
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Since AA f is the Cholesky parametrization of Pg , clearly K E Pg if and only if an > 0, i = 
1,2,3. If we can prove the following lemma, the condition M > 0 will yield the dual cone 

D* 

r Q' 

Lemma Al. The trace (K,D) is positive for all K E Pg \ {0} if and only if the matrix 
M of the quadratic form ( K , D) = a t Ma is positive definite 

Let us now prove the lemma. Clearly if M > 0 then ( K,D) = a 1 Ala > 0 for all a E R 5 
and in particular for all a with an > 0, * = 1,2,3. Conversely let a E R 5 . Then a can be 
written as 

a = (ei«.n, £ 2 ^ 22 , £ 3033 , ai 2 , 013 )* 
where e* is the sign of an, i = 1,2,3 and we have 

a 1 Ala = (a' 1 + af 2 + a i 3 )dn + (a 22 + 013)^22 + a 33^33 + 2 (e 2 <i 220 i 2 + 033)^12 + 463013033 ^ 13 . 
But this is also equal to a*Ala where 

h* = (|on|, |a 2 2|, I033I, 62012, 63013) 

which is in Pg. Therefore ( K,D) > 0 for all K E Pg if and only if M is positive definite 
which translates immediately into the conditions defining Pg in Proposition 14.31 


Proof of Theorem 14.41 For the proof of the theorem, it will be convenient to adopt a 
slightly different form of the parametrization of the Cholesky decomposition of K = AAL in 
Pg. Let 


Ai 



if i = j, 
if i < j. 


so that 


(AA T )ij 


an + a ii 

l>i 

Q'ij V Q'jj Q'il&jl 

l>max(i,j) 


if i = j, 
if i < j. 


Ecjuating each entry kij of K to the corresponding entry of ylyl T with the constraint that 


k\2 — —\Ja 22 a 12 + ai3«23) 


k 13 = k 23 shows that 
hi = an + af 2 + a? 3 , 
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&22 — «22 + ®23) k 13 — ~ V a 33 a 13) 

&33 = 033 , &23 = — v / “33 a 23- 

In particular, we find that since 033 > 0 , 013 = 023 and fci2 = —y/0,22(112 + af 3 . The 
Jacobian of the transformation from K to A is 


k\2 

0 

— \/ a 22 


kn 

an l 1 

G 12 * 

J = «13 * 

022 
«33 \ 


It is easy to see |J| = |dza<?(J)| = 0^22 ^33 ■ 


ki3 
0 

0 0 

~i/«33 2ai3 

* 1 

* * 


k22 ^33 

0 0 \ 


0 
0 
0 
1 


We now have all the ingredients necessary to calculate the normalizing constant Ig(5 : D ). 
We have \K\ = 011022033 and 

(K,D) = dnkn + ^ 22^22 + <^ 33^33 + 2 ^ 12^12 + 2dizki3 + 2 ^ 23^23 

= dii(an + Oi2 + of 3) + G?22 (022 + O23) + £^33033 

+2di2(-v / 0220i2 + 013023) + 2di 3 (—013-^/033) + 2 d 2 3 (— 023 v / 033 )- 

and so the normalizing constant is 

r 5-2 s-i s -1 j ^ 1 1 

Ig(5,D) = / a,f a 2 | a 33 2 exp(--dnan - -^22022 - 7:^33033 - -dna^ 

Ja z z z z 

~ 2^ U ^ 22 ^12) a 13 + ^12\/®22 a 12 + 2dl30i 3v /033)(i^4. 

where an > 0; E R,i < j; and cL4 denotes the product of all differentials. The integral 
with respect to an is a gamma integral with 

til 1 s S -l 

/ a n 2 exp(--dnaii)daii = 22 T(-)d n 2 . 

Jo 2 2 

The integral with respect to 012 and 013 are Gaussian integrals with 

/ exp(-idno 32 + di2v / ®22 a i2)dai2 = exp( d ^ 22 ), 

J -00 2 v 011 2an 
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and 


/ oo ^ 

exp(—-(dn + d 22 + 2di2)ai3 + 2di3^/a33ai3)dai3 = 

-OO ^ 

Therefore 

Ig(5,D) = r( —)2^d 11 2 27r(dn + d 22 + 2di2) 2 


+2T 


; P(- 


2df 3 a 3 3 


V dn + d 22 + 2 di 2 dn 4 d22 + 2 di 2 


r°° s-i s-i 4 d ^2 1 2df 3 

J a 22 a 33 exp ^~2 d22+ 2^^ a22 + ^~2 d33+ d^Td^T2d^^ a33 ^ da22da33 


— r( — )2 2 d 11 2 27r(dn + d22 + 2di2) 2 


p/ ^ + l w 2c hi + 2(dn + d 2 2 + 2di 2 ) ,i+l 

2 dnd22 — d \ 2 2 d33(dn + d 2 2 + 2di 2 ) — 4dj 3 

$ 0^+1 35+4 5 0 5+1 

= r(— )T (—-— )n 2 2 (dn + d>22 + 2^12) 2 [^33(^11 + d22 + 2d\2) — 4 g?i3 ] 2 

x (^11^22 — ^12) 2 • 


Proof of Proposition 14.41 

We proceed as in the proof of Proposition 14.31 That is, we let K = AA* be the Cholesky 
decomposition of K with + upper triangular. Equating the entries of K and AA t yields 


«23 — O43, «24 — 014, a 44 — 033 

with then 

2 2 2 2 2 2 
fell = Oil + a 12 + °12 + O44, k \2 = 012 022 + O43 + 0 14 , fcl3 = 013033. k \4 = 014033 

&22 = O22 + O43 + a 44 , k‘ 2 :i = 013033, fei4 = 014033 

&33 = O33 &34 = 0 

/C44 —— (Z33 

Then, ordering (K,D) as a polynomial in a^-, we see that 

(K,D) = dnah + < 3(22022 + 2 ^ 330^3 + dnaf 2 + 2^42022012 + 043(^11 + 2di2 + ^ 22 ) 

+4^13013033 + af 4 (dn + 2di2 + ^22) + 4^14014033 
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is a quadratic form and the matrix of this quadratic form is 


/ 

dn 

0 

0 

0 

0 

0 

\ 


0 

d -22 

di2 

0 

0 

0 



0 

di 2 

dn 

0 

0 

0 



0 

0 

0 

2 d 33 

2 <Ll 3 

2 di 4 



0 

0 

0 

2^43 

dn + 2g? 4 2 + d 22 

0 


V 

0 

0 

0 

2 c ?44 

0 

dn + 2 d 4 2 + d 22 / 


With exactly the same argument as in Proposition 14.31 we can show that (K, D) > 0 for all 
K E Pg if and only if M > 0 , i.e. D satisfies the conditions of Proposition 14.41 


Proof of Theorem 14.51 As in the proof of Theorem 14.41 it will be convenient to adopt a 
slightly different parametrization of the Cholesky decomposition of K. Let 


Aij — 


if i = j, 
if i < j. 


so that the entries of AAd are given by 
(AA T ) i;j = 


an + E a 2 


il 


—a. 


l>i 

V V a P ~*~ 


a u a ji 

z>max(ij) 


if i = j, 
if i < j. 


fell = an + a? 2 + a? 3 + a? 


Equating each entry kjj of K to the corresponding entry of AA T , we hnd that 

ki 2 = —\Ja22a\2 + 043023 + ai 4 a 2 4 , 

&14 = — i/a44ai4, 

^’23 = — v / °33 a 23 + 024034, 

&33 = 033 + a§ 4 , 

/C44 = 044. 

This shows that 044 > 0 and 034 = 0 . Since 033 > 0 and £43 = fc 2 3, then 043 = a 2 3. Since 
a 4 4 > 0 and /C44 = fc 2 4, then 044 = a 24 . Since /C34 = 0 , then 033 = 044. Therefore, we obtain 
that 

fell = 044 + a \ 2 + af 3 + af 4 , /c 12 = -^/di^a 12 + af 3 + af 4 , 

&13 = ^’23 = — \/ a 33 a 13) = fc 2 4 = ”1/033044, 


Hj 

*12 w a 13 w “14) 
&13 = ”i/ a 33 a 13 + 044034, 
&22 = 0 22 + 023 A a 24) 

&24 = —1/044 o 24 
&34 = —V / “ 44 Q ' 34 , 
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&22 — a 22 + ®13 + a 14 > 


&33 = &44 = O33. 


The Jacobian of the transformation from K to A is 


J = 


an 

012 

013 

014 

022 

033 


hi 

t 1 


k\2 

0 

\Z ®22 

* 

* 

* 

* 


h 3 

0 

0 

-■/033 

* 

* 

* 


k u 
0 

0 0 
0 2ai3 

— 1/033 2 ai 4 


&22 &33 

0 0 \ 
0 
0 
0 
0 

1 ) 


1 /2 

It is easy to see \J\ = \diag{J)\ = a 22 033. We now have all the ingredients necessary to 
calculate the normalizing constant /g 2 (< 5 , D). Through the change of variables, K = AA T . 
Then \K\ = 041022033, 

( K , D) = dllfcn -f- 0(22^22 + ^ 33^33 + d 44&44 + 2 dl 2&12 + 2 dl 3 &l 3 + 2 ^ 14^14 + 2 (^ 23^23 + 2 (^ 24^24 + 2 (^ 34^34 

= dn(an + aj2 + cti3 + 014) + ^22(022 + <Ji3 + 044) + 2^33033 

+2di2(—\/o 220 i 2 + 013 + a i 4 ) + 4 di 3 (—ai 3 \/a 33 ) + 4 di 4 (—ai 4 \/ 033). 


and so the integral equals 

-2 .5-1 




L 


L 11 <x 22 ^33 


1 1 2 1 2 
exp{—-dnan — —dna 12 + di2-/o22ffli2 — — (dn + d22 + 2di2)a 13 


+2di3ai3-/<233 — — (dll + d22 + 2dl2)ai4 + 2 dl 4 ai 4\/ a 33 — — d 2 20,22 — d33033}dT. 


where an > 0 ; € R,i < j\ and dA denotes the product of all differentials. The integral 

with respect to an is gamma integrals, then 

roo s _2 1 a 5 _s 

J a n 2 exp(— -diian)dan = 2 2 T(-)d 11 2 . 


The integral with respect to 042, 043 and 044 are normal integrals, then 


f 


exp(—7-^11042 + di 2 y /a^ai 2 )dai 2 = exp/^ 12022 ' 


y/dn 


2dn 


/oo 2 

exp(—-(dn + d22 + 2di2)a? 3 + 2di3-/a33ai3)dai3 

-OO ^ 


V2T 


exp(- 


2d 13 Cl33 


\/dn + *2 + 2*2 *1 H - *2 + 2*2 
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and 


F 


exp(-i(dn + d 2 2 + 2di 2 )o? 4 + 2di 4 v / aiJai4)dai 4 = v/ ~E ^^== exp( 

> ^ V “11 + “22 + 2di2 


2df 4 “33 


dll + d 2 2 + 2dl2 J 


Therefore, the integral becomes 

Igi (d, H) = r(-) 2 2 d n 2 (2^)5(du + d 2 2 + 2 di 2 ) 


L 


o° 5-1 5-1 

a 22 2 °33 2 exp{(--d 2 2 + 


2 0/7 2 _i_ 0 /7 2 

12 )“22 + ( — ^33 + —-—)“33}d“22d“33 


2dll 


dn + d 2 2 + 2di 2 ' 


5 a +3 _i±i 3 

= r(-)2 2 d 14 2 7r 2 (dn + d 2 2 + 2di 2 ) 

r ^+1 , 2 dn .^±1 . , ._ dn + d 2 2 + 2di2 _,5 

2 dnd 22 — dj 2 d 33 (dn + d 22 + 2 di 2 ) — 2 (d 43 + df 4 ) 

= r(|)r(^i)r( 5 ) 7 ri 2 ,s+2 (dii + d 22 + 2 di 2 )^ 1 [d 33 (dii + d 22 + 2 d i2 ) - 2 (di 3 + d ? 4 )]- 5 

(dnd 22 — d 2 2 ) 2 . 


Appendix 2 


Numerical values for D,E(K) and K and plots for Section 5 

We give here the matrices D, E{K) and K as well as the traceplot and autocorrelation 
plot of log(|if|) for all graphs in Fig. [T] Here K has been computed with 5000 iterations 
after a 1000 iterations burn in and averaged over 100 simulations 


Graph in Fig. [l](a) 
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Traceplot of Logdet 


Autocorrelation of Logdet 



(a) Traceplot (b) ACF plot 

Figure 3: (a) Traceplot of log(|i7|) v.s. the number of iterations, (b) Autocorrelation plot 
of log(|AT|) for Fig. Oja). 


1.1294 

0 

0 

-0.0129 

0 

0 

0 

0 

0.5915 

0 

-0.0129 

0 

0 

0 

0 

0 

0.2578 

-0.0129 

0 

0 

0 

-0.0129 

-0.0129 

-0.0129 

0.0767 

-0.0129 

0 

0 

0 

0 

0 

-0.0129 

0.2589 

-0.0129 

-0.0129 

0 

0 

0 

0 

-0.0129 

0.3699 

0 

0 

0 

0 

0 

-0.0129 

0 

0.2817 


1.1274 

0 

0 

-0.0127 

0 

0 

0 

\ 

0 

0.5961 

0 

-0.0127 

0 

0 

0 


0 

0 

0.2563 

-0.0127 

0 

0 

0 


-0.0127 

-0.0127 

-0.0127 

0.0767 

-0.0127 

0 

0 


0 

0 

0 

-0.0127 

0.2594 

-0.0127 

-0.0127 


0 

0 

0 

0 

-0.0127 

0.3708 

0 


0 

0 

0 

0 

-0.0127 

0 

0.2818 

/ 


36 


























Graph in Fig. [0(b) 
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Traceplot of Logdet 


Autocorrelation of Logdet 
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(b) ACF plot 

Figure 4: Traceplot and Autocorrelation plot of log(|iT|) for Graph in Fig. Htb). 
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Graph in Fig. [0(c) 
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(a) Traceplot (b) ACF plot 

Figure 5: Traceplot and Autocorrelation plot of log(|-fC|) for Graph in Fig. DJc). 
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Figure 6: Traceplot and Autocorrelation plot of log(|iT|) for Graph in Fig. [ijd) 
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Figure 7: Traceplot and Autocorrelation plot of log(|iC|) for Graph in Fig. [He) 
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Appendix 3 


Estimates and batch standard errors for entries of K for the models of 
Section 6 

The estimates and batch standard errors are given below for the entries of K listed in 
lexicographic order. 

Table 4: The average estimates for entries of K for Fig. [2] (b) when p = 20 


0.1072 

0.0100 

0.0096 

0.0322 

0.0109 

0.0093 

0.0102 

0.0105 

0.0101 

0.0103 

0.0099 

0.0106 

0.0100 

0.0099 

0.0109 

0.0104 

0.0111 

0.0116 

0.0104 

0.0100 

0.0113 

0.0115 







Table 5: The batch standard errors for Fig. [2] (b) when p = 20 


0.0004 

0.0005 

0.0005 

0.0001 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0005 

0.0004 

0.0005 

0.0004 

0.0005 

0.0004 

0.0005 

0.0005 

0.0005 

0.0005 

0.0004 

0.0005 




Table 6: The average estimates for entries of K for Fig. [2] (b) when p = 30 


0.1217 

0.0109 

0.0126 

0.0366 

0.0109 

0.0118 

0.0120 

0.0120 

0.0115 

0.0122 

0.0108 

0.0121 

0.0113 

0.0119 

0.0125 

0.0114 

0.0120 

0.0112 

0.0119 

0.0131 

0.0115 

0.0125 

0.0116 

0.0132 

0.0110 

0.0119 

0.0119 

0.0107 

0.0129 

0.0119 

0.0124 

0.0119 










Table 7: The 

batch standard errors for 

Fig. 0 (b) 

when p 

= 30 



0.0008 

0.0006 

0.0006 

0.0003 

0.0006 

0.0006 

0.0006 

0.0006 



0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 



0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 

0.0006 



0.0006 

0.0005 

0.0006 

0.0005 

0.0006 

0.0006 

0.0005 

0.0005 
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Table 8: The average estimates for entries of K for Fig. [2] (c) when p = 20 


0.1102 

0.0106 

0.0104 

0.0347 

0.1135 

0.0329 

0.1104 

0.0335 

0.1113 

0.0332 

0.1103 

0.0326 

0.1157 

0.0330 

0.1082 

0.0333 

0.1083 

0.0318 

0.1096 

0.0326 

0.1059 

0.0311 










Table 9: The 

batch standard < 

srrors for Fig. [2] (c) 

i when p ■ 

= 20 



0.0011 

0.0002 

0.0002 

0.0004 

0.0012 

0.0003 

0.0011 

0.0004 



0.0012 

0.0004 

0.0012 

0.0003 

0.0013 

0.0003 

0.0012 

0.0003 



0.0012 

0.0004 

0.0011 

0.0003 

0.0012 

0.0003 




Table 10: 

The average estimates for 

entries of 

I\ for Fig. [2] (c) when p = 

30 

0.1295 

0.0117 

0.0111 

0.0384 

0.1253 

0.0386 

0.1266 

0.0376 

0.1248 

0.0357 

0.1214 

0.0358 

0.1209 

0.0357 

0.1181 

0.0358 

0.1161 

0.0349 

0.1126 

0.0345 

0.1123 

0.0339 

0.1126 

0.0338 

0.1136 

0.0330 

0.1143 

0.0323 

0.1083 

0.0324 

0.1077 

0.0318 










Table 11: The 

batch standard 

errors for 

Fig. [2] (c) when p 

= 30 



0.0013 

0.0002 

0.0002 

0.0003 

0.0012 

0.0004 

0.0011 

0.0004 



0.0011 

0.0003 

0.0011 

0.0003 

0.0013 

0.0004 

0.0010 

0.0004 



0.0011 

0.0003 

0.0010 

0.0003 

0.0010 

0.0003 

0.0011 

0.0003 



0.0010 

0.0003 

0.0012 

0.0003 

0.0010 

0.0003 

0.0010 

0.0003 
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